load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
;======================================input===========================
  fread  = addfile("data/obs/precip/GPCC/full_data_monthly_v2022_05.nc","r")
  YYYYMM = cd_calendar(fread->time,1)
  tStrt  = ind(YYYYMM.eq.195101)
  tLast  = ind(YYYYMM.eq.202012)
  pre    = fread->pre(tStrt:tLast,{16:42},{65:105})

;=====================================calculation==============================
  pre!1 = "lat"
  pre!2 = "lon"

  pr     = month_to_season(pre,"JJA")
  pr     = (pre(5::12,:,:)+pre(6::12,:,:)+pre(7::12,:,:)+pre(8::12,:,:))/4.; summer(JJAS) mean

  pr     = runave_n_Wrap(pr,9,0,0)
  pr0    = pr ; original precipitation 
  clm    = conform(pr,dim_avg_n_Wrap(pr,0),(/1,2/))
  std    = conform(pr,dim_stddev_n_Wrap(pr,0),(/1,2/))
  pr     = where(std.eq.0,pr@_FillValue,pr)
  clm    = where(std.eq.0,clm@_FillValue,clm)
  std    = where(std.eq.0,std@_FillValue,std)
  pr     = (pr-clm)/std ; normalized precipitation anomalies

  rad      = 4.*atan(1.)/180.
  clat     = pr&lat
  clat     = sqrt(cos(rad*clat))
  wx       = pr
  wx       = pr*doubletofloat(conform(pr,clat,1))

  neof        = 2 
  optEOF      = True
  optEOF@jopt = 0
  optETS      = False

  x           = wx(lat|:,lon|:,time|:)
  eof         = eofunc_Wrap(x,neof,optEOF)
  eof_ts      = eofunc_ts_Wrap(x,eof,optETS)
  eof_ts      = dim_standardize_n(eof_ts,0,1); the corresponding time series of the first two leading modes
  eof_regress = regCoef(eof_ts,pr0(lat|:,lon|:,time|:)); related precipitation pattern
  copy_VarCoords(eof(:,0,0),eof_regress(:,0,0))
  copy_VarCoords(pr(0,:,:),eof_regress(0,:,:))
  ns          = dimsizes(pr)
  prinfo      = True
  sig         = eofunc_north(eof@pcvar,ns(0),prinfo); test

;============================output====================

  foutname = "data/pr_eof_1951-2020_JJAS_GPCC.nc"
  system("rm -f "+foutname)
  fout = addfile(foutname,"c")
  fout->eof_regress=eof_regress ; precipitation pattern
  fout->eof=eof ; with explained variances
  fout->ts=eof_ts; time series

end 


